function[y] = f1(p_star)
global gamma p1 rho1 c1 p2 rho2 c2 u1 u2;

if(p_star >= p1)
   y = (p_star - p1)/rho1/c1/( (gamma+1)/2/gamma *(p_star/p1) + (gamma-1)/2/gamma)^0.5;
end
if(p_star <p1)
    y = 2*c1/(gamma-1)*((p_star/p1)^((gamma-1)/2/gamma) - 1);
end

end